ATM 623: Climate Modeling¶
Brian E. J. Rose, University at Albany
Lecture 9: Who needs spectral bands? We do. Some baby steps...¶
About these notes:¶
This document uses the interactive IPython notebook format (now also called Jupyter). The notes can be accessed in several different ways:
- The interactive notebooks are hosted on
githubat https://github.com/brian-rose/ClimateModeling_courseware - The latest versions can be viewed as static web pages rendered on nbviewer
- A complete snapshot of the notes as of May 2015 (end of spring semester) are available on Brian's website.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
1. What if CO$_2$ actually behaved like a Grey Gas?¶
Suppose that CO$_2$ actually behaved as a grey gas. In other words, no spectral dependence in absorptivity.
If we then double the CO2 concentration in the atmosphere, we double the number of absorbers. This should imply that we also double the absorption cross-section:
$$ \kappa^\prime = 2 ~ \kappa $$which would imply that we double the optical thickness of every layer:
$$ \Delta \tau^\prime = 2 \left( -\frac{\kappa}{g} \Delta p \right) = 2 ~ \Delta \tau$$and since the absorptivity / emissivity of each layer is
$$ \epsilon = 1 - \exp\big( - \Delta \tau \big) $$the modified absorptivity is
$$ \epsilon^\prime = 1 - \exp\big( - 2\Delta \tau \big) = 1 - \left( \exp\big( - \Delta \tau \big)\right)^2 = 1 - (1-\epsilon)^2 $$or simply $$ \epsilon^\prime = 2 \epsilon - \epsilon^2 $$
(Note that $\epsilon^\prime = 2 \epsilon$ for very thin layers, for which $\epsilon$ is small).
What does our 2-layer analytical model then say about the radiative forcing?¶
Recall that we tuned the two-layer model with
$$ \epsilon = 0.58 $$to get the observed OLR with observed temperatures.
If CO2 behaved like a grey gas, doubling it would then change the absorptivity from 0.58 to 0.82 (applying the above formula).
That's a 41% increase in the absorptivity of each layer!
Back in Lecture 6 we worked out that the radiative forcing in this model (with the observed lapse rate) is about +2.2 W m$^{-2}$ for a 1% increase in $\epsilon$.
This means that our hypothetical doubling of "grey CO$_2$" should yield a radiative forcing of 90.2 W m$^{-2}$.
This is an absolutely enormous number. Assuming a net climate feedback of -1.3 W m$^{-2}$ K$^{-1}$ (consistent with the AR5 ensemble) would then give us an equilibrium climate sensitivity of 70 K.
Conclusions:¶
- If CO2 did behave like a grey gas, we would be toast.
- The Grey Gas model is insufficient for understanding radiative forcing and feedback.
2. Another look at observed spectra¶
It's time to move away from the Grey Gas approximation and look more carefully at the actual observed spectra of solar and terrestrial radiation.
Observed solar spectra¶
The following figure shows observed spectra of solar radiation at TOA and at the surface, along with the theoretical Planck function for a blackbody at 5525 K.
from IPython.display import Image
Image('../images/Solar_spectrum.png')
This figure shows the solar radiation spectrum for direct light at both the top of the Earth's atmosphere and at sea level. The sun produces light with a distribution similar to what would be expected from a 5525 K (5250 °C) blackbody, which is approximately the sun's surface temperature. As light passes through the atmosphere, some is absorbed by gases with specific absorption bands. Additional light is redistributed by Raleigh scattering, which is responsible for the atmosphere's blue color. These curves are based on the American Society for Testing and Materials (ASTM) Terrestrial Reference Spectra, which are standards adopted by the photovoltaics industry to ensure consistent test conditions and are similar to the light that could be expected in North America. Regions for ultraviolet, visible and infrared light are indicated.
Source: http://commons.wikimedia.org/wiki/File:Solar_spectrum_en.svg
The figure shows that that the incident beam at TOA has the shape of a blackbody radiator.
By the time the beam arrives at the surface, it is strongly depleted at specific wavelengths.
Absorption by O$_3$ (ozone) depletes almost the entire ultraviolet spectrum.
Weaker absorption features, mostly due to H$_2$O, deplete some parts of the near-infrared.
Note that the depletion in the visible band is mostly due to scattering, which depletes the direct beam but contributes diffuse radiation (so we can still see when it's cloudy!)
Observed terrestrial spectra¶
This figure shows the Planck function for Earth's surface temperature compared with the spectrum observed from space.
Image('../images/Terrestrial_spectrum.png')
Careful: I'm pretty sure what is plotted here is not the total observed spectrum, but rather the part of the emissions from the surface that actual make it out to space.
As we now, the terrestrial beam from the surface is depleted by absorption by many greenhouse gases, but principally CO$_2$ and H$_2$O.
However there is a spectral band centered on 10 $\mu$m in which the greenhouse effect is very weak. This is the so-called window region in the spectrum.
Since absorption is so strong across most of the rest of the infrared spectrum, this window region is a key determinant of the overall greenhouse effect.
One very big shortcoming of the Grey Gas model: it ignores the window region¶
We would therefore like to start using a model that includes enough spectral information that it represents
- the mostly strong CO2 absorption outside the window region
- the weak absorption inside the window region
3. Water vapor changes under global warming¶
Another big shortcoming of the Grey Gas model is that it cannot represent the water vapor feedback.
We have seen above that H$_2$O is an important absorber in both longwave and shortwave spectra.
We also know that the water vapor load in the atmosphere increases as the climate warms. The primary reason is that the saturation vapor pressure increases strongly with temperature.
Evidence from CESM simulations¶
Let's take at changes in the mean water vapor fields in the CESM model after a doubling of CO$_2$
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
# Open handles to the data files
datapath = "http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/Brian+Rose/CESM+runs/"
endstr = "/entry.das"
ctrl = nc.Dataset( datapath + 'som_control/som_control.cam.h0.clim.nc' + endstr )
co2 = nc.Dataset( datapath + 'som_2xCO2/som_2xCO2.cam.h0.clim.nc' + endstr )
# Get grid dimensions
lat = ctrl.variables['lat'][:]
lon = ctrl.variables['lon'][:]
lev = ctrl.variables['lev'][:]
# Get temperature, specific humidity and relative humidity
print ctrl.variables['T']
print ctrl.variables['Q']
print ctrl.variables['RELHUM']
t_ctrl = ctrl.variables['T'][:]
q_ctrl = ctrl.variables['Q'][:]
rh_ctrl = ctrl.variables['RELHUM'][:]
t_co2 = co2.variables['T'][:]
q_co2 = co2.variables['Q'][:]
rh_co2 = co2.variables['RELHUM'][:]
# Average in time and longitude
t_ctrl_zon = np.mean(t_ctrl, axis=(0,3))
t_co2_zon = np.mean(t_co2, axis=(0,3))
q_ctrl_zon = np.mean(q_ctrl, axis=(0,3))
q_co2_zon = np.mean(q_co2, axis=(0,3))
rh_ctrl_zon = np.mean(rh_ctrl, axis=(0,3))
rh_co2_zon = np.mean(rh_co2, axis=(0,3))
fig = plt.figure(figsize=(14,4))
axlist = []
ax = fig.add_subplot(1,3,1)
CS = ax.contourf(lat, lev, t_co2_zon - t_ctrl_zon,
levels=np.arange(-11,12,1), cmap=plt.cm.seismic)
ax.set_title('Temperature (K)')
fig.colorbar(CS)
axlist.append(ax)
ax = fig.add_subplot(1,3,2)
CS = ax.contourf(lat, lev, (q_co2_zon - q_ctrl_zon)*1000,
levels=np.arange(-3,3.25,0.25), cmap=plt.cm.seismic)
ax.set_title('Specific humidity (g/kg)')
fig.colorbar(CS)
axlist.append(ax)
ax = fig.add_subplot(1,3,3)
CS = ax.contourf(lat, lev, (rh_co2_zon - rh_ctrl_zon),
levels=np.arange(-11,12,1), cmap=plt.cm.seismic)
ax.set_title('Relative humidity (%)')
fig.colorbar(CS)
axlist.append(ax)
for ax in axlist:
ax.invert_yaxis()
ax.set_xticks([-90, -60, -30, 0, 30, 60, 90]);
ax.set_xlabel('Latitude')
ax.set_ylabel('Pressure')
plt.show()
A few things to see here:
- The warming is largest at
- the polar surface
- the tropical upper troposphere
- The specific humidity anomaly (i.e. the moistening) is
- largest at the tropical surface
- decreases upward and poleward
- In fact the specific humidity anomaly has roughly the same shape of the specific humidity field itself -- it is largest where the temperature is highest. This is a consequence of the Clausius-Clapeyron relation.
- The relative humidity anomaly is
- overall rather small (just a few percent)
- Largest in places cold places where the specific humidity is very small.
The smallness of the relative humidity change is a rather remarkable result.
This is not something we can derive from first principles. It is an emergent property of the GCMs. However it is a very robust feature of global warming simulations.
4. A simple water vapor parameterization¶
A credible climate model needs a water vapor feedback¶
If relative humidity is nearly constant under global warming, and water vapor is a greenhouse gas, this implies a positive feedback that will amplify the warming for a given radiative forcing.
Thus far our simple models have ignored this process, and we have not been able to use them to assess the climate sensitivity.
To proceed towards more realistic models, we have two options:
- Simulate all the evaporation, condensation and transport processes that determine the time-mean water vapor field (as is done in the CESM).
- Parameterize the dependence of water vapor on temperature by insisting that relative humidity stays constant as the climate changes.
We will now explore this second option, so that we can continue to think of the global energy budget under climate change as a process occurring in a single column.
Manabe's constant relative humidity parameterization¶
We are going to adopt a parameterization first used in a very famous paper:
Manabe, S. and Wetherald, R. T. (1967). Thermal equilibrium of the atmosphere with a given distribution of relative humidity. J. Atmos. Sci., 24(3):241–259.
This paper was the first to give a really credible calculation of climate sensitivity to a doubling of CO2 by accounting for the known spectral properties of CO2 and H2O absorption, as well as the water vapor feedback!
The parameterization is very simple:
We assume that the relative humidity $r$ is a linear function of pressure $p$:
$$ r = r_s \left( \frac{p/p_s - 0.02}{1 - 0.02} \right) $$where $p_s = 1000$ hPa is the surface pressure, and $r_s$ is a prescribed surface value of relative humidity. Manabe and Wetherald set $r_s = 0.77$, but we should consider this a tunable parameter in our parameterization.
Since this formula gives a negative number above 20 hPa, we also assume that the specific humidity has a minimum value of $0.005$ g/kg (a typical stratospheric value).
This formula is implemented in climlab.radiation.water_vapor.ManabeWaterVapor()
Using this parameterization, the surface and tropospheric specific humidity will always increase as the temperature increases.
Here is a brief introduction to the climlab.BandRCModel process.
This is a model that divides the spectrum into 7 distinct bands: three shortwave and four longwave.
As we will see, the process works much like the familiar climlab.RadiativeConvectiveModel.
About the spectra¶
The shortwave is divided into three channels:
- Channel 0 is the Hartley and Huggins band (extreme UV, 200 - 340 nm, 1% of total flux, strong ozone absorption)
- Channel 1 is Chappuis band (450 - 800 nm, 27% of total flux, moderate ozone absorption)
- Channel 2 is remaining radiation (72% of total flux, largely in the visible range, no ozone absorption)
The longwave is divided into four bands:
- Band 0 is the window region (between 8.5 and 11 $\mu$m), 17% of total flux.
- Band 1 is the CO2 absorption channel (the band of strong absorption by CO2 around 15 $\mu$m), 15% of total flux
- Band 2 is a weak water vapor absorption channel, 35% of total flux
- Band 3 is a strong water vapor absorption channel, 33% of total flux
The longwave decomposition is not as easily related to specific wavelengths, as in reality there is a lot of overlap between H$_2$O and CO$_2$ absorption features (as well as absorption by other greenhouse gases such as CH$_4$ and N$_2$O that we are not representing).
Example usage of the spectral model¶
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
from climlab import constants as const
First try a model with all default parameters. Usage is very similar to the familiar RadiativeConvectiveModel.
col1 = climlab.BandRCModel()
print col1
Check out the list of subprocesses.
We now have a process called H2O, in addition to things we've seen before.
This model keeps track of water vapor. We see the specific humidity in the list of state variables:
col1.state
The water vapor field is initialized to zero. The H2O process will set the specific humidity field at every timestep to a specified profile. More on that below. For now, let's compute a radiative equilibrium state.
col1.integrate_years(2)
# Check for energy balance
col1.diagnostics['ASR'] - col1.diagnostics['OLR']
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot( col1.Tatm, col1.lev, 'c-', label='default' )
ax.plot( col1.Ts, climlab.constants.ps, 'co', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_title('Temperature profiles', fontsize = 18)
ax.grid()
By default this model has convective adjustment. We can set the adjusted lapse rate by passing a parameter when we create the model.
The model currently has no ozone (so there is no stratosphere). Not very realistic!
More reasonable-looking troposphere, but still no stratosphere.
About the radiatively active gases¶
The Band model is aware of three different absorbing gases: O3 (ozone), CO2, and H2O (water vapor). The abundances of these gases are stored in a dictionary of arrays as follows:
col1.absorber_vmr
Ozone and CO2 are both specified in the model. The default, as you see above, is zero ozone, and constant (well-mixed) CO2 at a volume mixing ratio of 3.8E-4 or 380 ppm.
Water vapor is handled differently: it is determined by the model at each timestep. We make the following assumptions, following a classic paper on radiative-convective equilibrium by Manabe and Wetherald (J. Atmos. Sci. 1967):
- the relative humidity just above the surface is fixed at 77% (can be changed of course... see the parameter
col1.relative_humidity - water vapor drops off linearly with pressure
- there is a small specified amount of water vapor in the stratosphere.
Putting in some ozone¶
We need to provide some ozone data to the model in order to simulate a stratosphere. As we did with the original column model, we will use the ozone climatology data provided with the CESM model.
See here for more information, including some plots of the ozone data: http://www.atmos.albany.edu/facstaff/brose/classes/ENV480_Spring2014/styled-5/code-3/index.html
import netCDF4 as nc
datapath = "http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/Brian+Rose/CESM+runs/"
endstr = "/entry.das"
topo = nc.Dataset( datapath + 'som_input/USGS-gtopo30_1.9x2.5_remap_c050602.nc' + endstr )
ozone = nc.Dataset( datapath + 'som_input/ozone_1.9x2.5_L26_2000clim_c091112.nc' + endstr )
# Dimensions of the ozone file
lat = ozone.variables['lat'][:]
lon = ozone.variables['lon'][:]
lev = ozone.variables['lev'][:]
# Taking annual, zonal, and global averages of the ozone data
O3_zon = np.mean( ozone.variables['O3'],axis=(0,3) )
O3_global = np.sum( O3_zon * np.cos(np.deg2rad(lat)), axis=1 ) / sum( np.cos(np.deg2rad(lat) ) )
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot( O3_global*1E6, lev)
ax.invert_yaxis()
We are going to create another instance of the model, this time using the same vertical coordinates as the ozone data.
# Create the column with appropriate vertical coordinate, surface albedo and convective adjustment
col2 = climlab.BandRCModel(lev=lev)
print col2
# Set the ozone mixing ratio
# IMPORTANT: we need to flip the ozone array around because the vertical coordinate runs the wrong way
# (first element is top of atmosphere, whereas our model expects the first element to be just above the surface)
col2.absorber_vmr['O3'] = np.flipud(O3_global)
# Run the model out to equilibrium!
col2.integrate_years(2.)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot( col1.Tatm, np.log(col1.lev/1000), 'c-', label='RCE' )
ax.plot( col1.Ts, 0, 'co', markersize=16 )
ax.plot(col2.Tatm, np.log(col2.lev/1000), 'r-', label='RCE O3' )
ax.plot(col2.Ts, 0, 'ro', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('log(Pressure)', fontsize=16 )
ax.set_title('Temperature profiles', fontsize = 18)
ax.grid()
ax.legend()
Once we include ozone we get a well-defined stratosphere. We can also a slight cooling effect in the troposphere.
Things to consider / try:
- Here we used the global annual mean Q = 341.3 W m$^{-2}$. We might want to consider latitudinal or seasonal variations in Q.
- We also used the global annual mean ozone profile! Ozone varies tremendously in latitude and by season. That information is all contained in the ozone data file we opened above. We might explore the effects of those variations.
- We can calculate climate sensitivity in this model by doubling the CO2 concentration and re-running out to the new equilibrium. Does the amount of ozone affect the climate sensitivity? (example below)
- An important shortcoming of the model: there are no clouds! (that would be the next step in the hierarchy of column models)
- Clouds would act both in the shortwave (increasing the albedo, cooling the climate) and in the longwave (greenhouse effect, warming the climate). Which effect is stronger depends on the vertical structure of the clouds (high or low clouds) and their optical properties (e.g. thin cirrus clouds are nearly transparent to solar radiation but are good longwave absorbers).
col3 = climlab.process_like(col2)
print col3
# Let's double CO2.
col3.absorber_vmr['CO2'] *= 2.
col3.compute_diagnostics()
print 'The radiative forcing for doubling CO2 is %f W/m2.' % (col2.diagnostics['OLR'] - col3.diagnostics['OLR'])
col3.integrate_years(3)
col3.diagnostics['ASR'] - col3.diagnostics['OLR']
print 'The Equilibrium Climate Sensitivity is %f K.' % (col3.Ts - col2.Ts)
col4 = climlab.process_like(col1)
print col4
col4.absorber_vmr
col4.absorber_vmr['CO2'] *= 2.
col4.compute_diagnostics()
print 'The radiative forcing for doubling CO2 is %f W/m2.' % (col1.diagnostics['OLR'] - col4.diagnostics['OLR'])
col4.integrate_years(3.)
col4.diagnostics['ASR'] - col4.diagnostics['OLR']
print 'The Equilibrium Climate Sensitivity is %f K.' % (col4.Ts - col1.Ts)
Interesting that the model is MORE sensitive when ozone is set to zero.
Credits¶
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences, offered in Spring 2015.
%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information
%version_information numpy, climlab